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Abstract 

Consider the macroscale modelling of microscale spatiotemporal 
dynamics. Here we develop a new approach to ensure coarse scale 
discrete models preserve important self-adjoint properties of the fine 
scale dynamics. The first part explores the discretisation of microscale 
continuum dynamics. The second addresses how dynamics on a fine 
lattice are mapped to lattice a factor of two coarser (as in multigrids). 
Such mapping of discrete lattice dynamics may be iterated to empower 
us in future research to explore scale dependent emergent phenomena. 
The support of dynamical systems, centre manifold, theory ensures 
that the coarse scale modelling applies with a finite spectral gap, in a 
finite domain, and for all time. The accuracy of the models is limited 
by the asymptotic resolution of subgrid coarse scale processes, and is 
controlled by the level of truncation. As given examples demonstrate, 
the novel feature of the approach developed here is that it ensures the 
preservation of important conservation properties of the microscale 
dynamics. 
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1 Introduction 

Dynamical systems theory gives new assurances about the quality of finite 
difference and finite element models of nonlinear spatiotemporal systems. Ef- 
forts to construct approximations to the long-term, low- dimensional dynam- 
ics of dissipative partial differential equations (pdes), on its inertial mani- 
fold [IHl e.g.], have largely been based upon the global nonlinear Galerkin 
method [29], [TH [25], e.g.], and its variants [I9l[l5l e.g.]. In contrast, a 'holistic 
discretisation' [30] developed further here is based purely upon the local dy- 
namics on finite elements while maintaining, as do inertial manifolds, fidelity 
with the solutions of the original pde. 
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To generate a macroscale model of a dissipative pde, we divide space 
into finite elements, then specially crafted coupling conditions empower us 
to support a macroscale model with centre manifold theory |8l[22], e.g.]. Such 
macroscale models are used for computations and for analytic understanding 
of the dynamics. Crucially, the theory supports the existence, exponentially 
quick attractiveness, and approximate construction of a slow manifold of the 
system dynamics, both deterministic [501 e.g.] and stochastic Earlier 
research developed coupling conditions that both had centre manifold sup- 
port on finite elements and ensured consistency as the element size became 
small [331 [211 e.g.]. However, albeit effective for many systems, these coupling 
conditions fail to preserve at each level of approximation any self-adjoint sym- 
metry in the underlying spatial dynamics [HT] . The innovations introduced in 
Section [2] are interelement coupling conditions that not only engender centre 
manifold support. Section [3l and assure consistency for vanishing element 
size, but also preserve self-adjoint symmetry in each approximation. For just 
one example. Section |l]2l as deduced in the model fl2T|) . argues that a partic- 
ular one-dimensional, nonlinear, continuum diffusion is soundly mapped as 
follows to nonlinear dynamics on a macroscale grid (spacing h): 
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where the grid values Uj = u(Xj,t). Why this discretisation instead of 
others? Because the new coupling conditions as well as having dynamical 
systems support and classic consistency, also preserve in the model the self- 
adjoint symmetry of material conservation that is present in the original 

PDE. 

In many science and engineering applications it is essential to preserve the 
conservative form of the governing equations. I argue that interelement cou- 
pling rules akin to those of Section [2] will automatically preserve conservative 
forms. 

Most methods for modelling dynamics posit just two time scales: a fast 
and a slow scale [HI [131 [2H1 e.g.]. Indeed, Sections [2H1| implicitly separate 
the dynamics of dissipative pdes into the 'uninteresting' fast subgrid dynam- 
ics, and the relevant slow dynamics of macroscale evolution resolved by the 
discretisation. But many applications possess a wide variety of interesting 
space-time scales [SI [121 e.g.]. Recent research developed a methodology with 
rigorous support for changing the resolved spatial grid scale by just a factor 
of two [32] • Homogenisation, in Section [721 is one example: the derivation of 
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equation fl5Up recommends that the evolution of discrete diffusion on a grid, 
with spatially varying diffusivity, is mapped to a coarser grid as 



^ ^^l^{^MUj-2-(/Cj_i+/Cj+i]Uj + /CHi^+2}, (2) 

where the coarse grid index j = 21 and the coarser scale diffusivity /Cj ~ 
^(K2j-2 + K2j-i + K2j+i + K2j+2) is an average over the microgrid diffusivi- 
ties. The mapping of dynamics from a finer grid to a coarser grid, via finite 
elements formed from a small number of fine grid points, may then be iter- 
ated to generate a hierarchy of models across a wide range of spatial scales, 
with the theory of centre manifolds to support across the whole hierarchy. 
This approach promises to empower us with great flexibility in modelling 
complex dynamics over multiple scales. Sections [5H7] further develop this 
modelling transformation of discrete dynamics on grids by exploring cou- 
pling conditions which result in the coarse grid dynamics also preserving the 
self-adjoint symmetries of the fine grid dynamics. 

Most two scale modelling methods can also be applied over many scales. 
However, most established methods require each application to be based upon 
a large 'spectral gap': a parameter such as e measures the scale separation, 
and invoking "as e — > 0" provides the extreme scale separation. In contrast, 
multigrid iteration for solving linear equations transforms between length 
scales that are different by (usually) a factor of two [3 [32], [6], e.g.]. Anal- 
ogously, Sections [5H7| explore modelling dynamics on a hierarchy of length 
scales that differ by a factor of two and hence the 'spectral gap' is finite, not 
infinite as required by popular extant, non-multigrid, methods for modelling 
dynamics. Section [H] describes how to divide fine grid lattice dynamics into 
small finite elements, develops interelement coupling rules that preserve self- 
adjoint symmetries, and then establishes the centre manifold support for the 
resulting coarse grid models. Section [7] outlines three applications including 
the homogenisation ([2]). 

The methodology proposed here uses dynamical systems theory to sup- 
port and construct accurate coarse grid models of fine scale dynamics, both 
continuum and discrete. I expect that the basis developed here for deter- 
ministic dynamics in one spatial dimension can be extended to both higher 
dimensions and stochastic dynamics. For non-dissipative dynamics, although 
the formal construction of coarse models follows analogously, current theory 
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Figure 1: to discretise continuum dynamics, bottom, to the grid, top, rewrite 
the continuum dynamics of u(x, t) as the dynamics of Uj(x, t) on overlapping 
elements. This figure plots three consecutive elements (blue, black and ma- 
genta): the jth element stretches from Xj_i to Xj+i. Then analysis supports 
and generates rules for the evolution of grid values U^(t) = Uj(Xj,t). 



gives little support for the relevance of the resulting sub-centre slow manifold 
models |^ [2l e.g.]. This article is confined to one dimensional dissipative 
dynamics. 

I also conjecture an implication for the equation-free modelling methodol- 
ogy of Kevrekidis et al. [201 [23l [T71 HOl e.g.]. When coupling sparse patches of 
microsimulators, previous research established that an analogous dynamical 
systems approach provided coupling of patches with high order accuracy on 
the macroscale [371 EH] • For dynamical systems where we want to automat- 
ically preserve in the macroscale any microscale symmetries, the coupling 
condition ([5|), shown to be required here, suggests that each patch should 
have a point source at mid-patch in proportion to fluxes extrapolated from 
neighbouring patches. Further research will tell. 



2 Self-adjoint preserving coupling conditions 

The next three sections explore the dynamics of a continuum field u(x, t) 
in one spatial dimension. Figure [1] shows a part of an L-periodic domain, 
u(x + L, t) — u(x, t], which is divided into m overlapping elements Fj = {x | 
Xj-i < X < Xj+i) m grid points Xj. These grid points need not be equally 
spaced. Typically Uj(x) and Vj(x) denote fields in the jth element Ej. Define 
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the inner product 



) = ^ uv dx = ^ 
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Xj-1 J 


J 



(3) 



These integrals in the inner product are spht over the two halves of each 
element to emphasise that perhaps unexpected contributions come from the 
central grid point Xj. 

Theorem 1 (self-adjoint coupling) The operator L in the linear system 
of coupled PDEs 
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1 > 
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(4) 



is self-adjoint when coupled with the interelement coupling conditions 

fj(X+)ujJX+)-fj(XpUjJXp 
+ y [f j+i (X^Uj+i ,,(Xj ) - f j_i (Xj )uj_i ,,(Xj )] 
- y ' [f j (Xj+i )Ujx(Xj+i ) - f j (Xj_i )Uj JXj_i )] = , 
Uj(Xj±i) =y'uj(Xj) +yUj±i(Xj±i) anc? Uj(Xr) = Uj(Xp , 

in which subscript x denotes spatial differentiation. 



(5) 
(6) 



Usually I link the coupling parameters by y + y' — 1 ; however, this 
constraint is not necessary in this theorem. 

In many applications the coefficient functions fj and do not depend 
upon the element j. However, in nonlinear systems we may need the reas- 
surance of the theorem when applied with nonlinear f(u, Ux] whence fj = 
f(Uj,Ujx) will be element dependent. Any dependence upon time and other 
parameters are suppressed for clarity. 

Proof: Undergraduate algebra proves the theorem. Integration by parts 
gives 

{Cm, v) = I [f jUjxVj - f jUjVj J + [f jUjxVj - f jUjV^ J 



+ 



Uj£vj dx 



(which upon using becomes) 

= (u, Cv) + Y_ {fj(Xr)uj,(Xr)vi(Xj] - i^[X;]u^[X^]v^^[Xr) 
i 

- f j (Xj_i )uj.(Xj_i ) [y 'v^ (Xj ) + yvj_i (Xj_i ]] 
+ f 5 (Xj_i ] [y 'Uj (X^ ) + yuj_i (X^_i )] VjJX^_i ) 
+ f 5 (Xj+i )uj.(Xj+i ] [y 'vj (Xj ] + yvj+i (X^i )] 

- f j (Xj+i ) [y 'uj (Xj ) + yuj+i (X^i )] Vj^X^i ) 
-f5(X,+]u^JX,+]vj(X3)+fj(X+)uj(X^^,(X+)} 

(upon renumbering Uj±^{Xj±■\] by j' = j ± 1 becomes) 

= (u, Cv) + Y_ {uj(X^) [ - f^XrjvjJXr) + y'fj(Xj_i)vjx(X^-i] 
j 

+ yf j+i (Xj ]vj+i ,,(Xj ) - y 'f ^ (Xj+1 IvjJX^i ) 
-yfj_i(Xj)vj_i,,(X3)+f^(Xt]v^,(X+]] 
+ Vj (Xj] [f j (Xr)uj,(Xr) - y% (Xj_i )ui.(Xj_i ) 
- yf j+1 (Xj ]uhi ,x(Xj ] + y 'f j (X^+i )uj JX^+i ) 
+ yfj_i(Xj)u3-i,.(Xj)-fj(X,+)uj,,(X^)]} 
= (u,£v) + by©. 

Hence £ with couphng conditions (l5])-(l6]) is self-adjoint. This proof also ap- 
plies, with only minor modifications, to multi-component systems where the 
field Uj in each element is a vector and f j and are symmetric matrices. ^ 



Remark The virtue of such coupling conditions for overlapping domains 
is that when finding the adjoint, the algebra involves quantities that are 
already involved, namely the grid values. If one tries to find the adjoint 
for disjoint elements, then the field values at the element boundaries also 
are involved, thus leading to too many free variables after the integration 
by parts. The overlap of finite elements used herein are increasingly being 
used in multiscale modelling: examples include the 'border regions' of the 
heterogeneous multiscale method [131 e-g-]? the 'buffers' of the gap-tooth 
scheme |40l e.g.], and the overlapping domain decomposition that improves 
convergence in waveform relaxation of parabolic pdes jl6[ e.g.]. 
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3 Centre manifold theory supports discreti- 
sation 



Based upon the equilibria and spectra about equilibria, centre manifold the- 
ory rigorously supports the existence, relevance and construction of low di- 
mensional models of dynamical systems [Bl [221 e-g-]- This section explores 
the theoretical support for forming discretisations of the class of self-adjoint 
PDEs in the form 



9u 
9t 



3 

9x 



f(x,U,Ux) — 

ox 



agix,u,Ux 



(7) 



using the self-adjoint preserving coupling conditions ©-([G]) with y' — 1 — y . 

Embed the L-periodic dynamics of the pde ([7]) into the overlapping ele- 
ments of Figure [1] by introducing the subgrid field Uj(x, t) in each element Ej. 
Define these subgrid fields to satisfy the pde ([7]) in each element, namely 



9uj 



8 

dx 



f(X,Uj,Ujx)-g^ 



+ ag(x,Uj,Ujx) , j = 1,...,m, (8) 



. ,Ura, y, oc) there exist a sub- 
a = (y' = 1) and subgrid 



and coupled by the conditions ([5])-([6]). 

In the extended state space E = (ui, 
space Eo of equilibria with parameters y : 
fields constant in each element, Uj(x, t] = Uj ; that is, fields which are piece- 
wise constant over the domain are equilibria when y = a = . For each 
element define fj(x) = f(x,Uj,0) then the differential equation for pertur- 
bations u-[x,t], linearised about the piecewise constant equilibria, is simply 
the diffusion equation 



9u; 



a 

ax 



1 



1 



.,m, 



(9) 



with, as y = (y' = 1), the 'insulating' version of the coupling condi- 
tions ©-([n]), namely for j = 1 , . . . , m 



fj(Xt)u;jXt)-fj(Xr)u;jXp 

- f J (Xj+1 )u;jXj+i ) + f j (Xj_i )u;jXj_i ; 
u;(Xj±o=u;(xr]=u;(x,+). 



0, 



(10) 

(11) 
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The centre manifold support for discrete models of the nonlinear sys- 
tem ([7]) with coupled elements rests upon the eigenstructure of the linearised 
system fl9l)- f|TT]) . As for all self-adjoint systems, following the proof of orthog- 
onality of eigenmodes, Section [3^ I proceed to prove that the eigenvalues are 
real and they are all non-positive. The proofs are elementary undergraduate 
algebra. Then Section 13.31 uses these to prove the existence and relevance of 
a slow manifold, 'holistic' discretisation for the nonlinear system ([7]). 

3.1 Homogeneous spectrum on an equi-spaced grid 

As one important example, this subsection considers the special case where 
the grid is equispaced and the 'diffusivity' f only depends upon space through 
its dependence upon the field u, that is, f = f(u, Ux). Then the linearised 
dynamics about the piecewise constant solutions are described by the pde ([9]) 
with coefficients fj which are constant on each element. Assume all fj are 
bounded above zero. This subsection then finds the spectrum and eigenmodes 
of the linearised pde 0. 

Seek eigenvalues A in the jth element such that fjU^^ = Au' (dropping 
subscript j for simplicity). As the eigenvalues must be real (see Theorem |3]), 
set A = — fjk^ for some k(> 0) to be determined. For an equi-spaced grid, 
AXj = h, , solutions must be of the form 

u' = Acosk(x-Xj] -hBsink(x-Xj) Csink|x - Xj| , (12) 

upon using the continuity at x = Xj . Then the insulating boundary condi- 
tions (dH) of u'(Xj±i) = u'(Xj) require A(coskh - 1) + (C ± B)sinkh = 0. 
The difference of these two conditions is 2B sin kh, = and hence either B = 
or kh, — nn . Explore both in turn. 

kh = TLTt Then coskh, = (—1)"' and the two conditions reduce to simply 
Afl-l)"^ - 1] = for which two cases arise depending upon even or 
odd values for n. 

even n Then A, B and C are unrestrained except possibly by the 
derivative condition in ffTOl) . However, all modes independently 
satisfy the derivative condition ffTOl) . Hence the three orthogonal 
modes u' = cos k(x — Xj ) , u' = sin k(x — Xj ] and u' — sin k|x — Xj | 
form a basis for the three dimensional eigenspace corresponding 
to eigenvalue A — — f jk^ . 
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odd n Then A = . Furthermore, sink(x — Xj) satisfies the derivative 
condition in fllOp . but sink|x — Xj| does not, C = . Hence u' = 
sink(x — X^) is the only eigenmode. 

One may hke to view the modes sinkjx — Xj| for even n with 
wavenumber k — nn/h as serving in place of the usual Fourier 
modes cosk(x — Xj) for odd n. 

kh ^ TLTT In this case, necessarily B = 0; then the condition A(coskh,— 1 ) + 
C sin kh = and the derivative condition flTU]) become 

cos kh. — 1 sin kh 
— sin kh cos kh — 1 

The determinant (coskh — 1)^ + sin^kh = 2(1 — coskh) is non-zero 
except for the outlawed cases kh = 2nn . Hence, A = C = also, and 
so there are no further eigenmodes. 

The spectrum of the linearised dynamics flU])- ffTT]) (and multiplicity) on 
the jth element is thus 

{0, -f jTrVH^ -4f jTrVH^ (triple) , -9fj7rVH^ -1 6i^n^/h^ (triple), . . .} . (13) 

This set, with one zero eigenvalue and the rest negative, provided 'diffusivity' 
required for the support of centre manifold theory. 

3.2 General spectral properties 

This subsection proves, using elementary undergraduate algebra, that crucial 
properties of the spectrum f|T3|) hold for the more general system fl9l)- f|TT]) 
for possibly non-uniform elements. Alternatively, one may view the first 
two of these properties as a consequence of the Self-adjoint Theorem [1], and 
the third as a natural consequence of the dissipation of diffusion. A reader 
comfortable with such a view may proceed direct to the next Section 13.31 

Throughout this subsection we only address a generic jth element in iso- 
lation, the interelement coupling parameter y — . This isolation is to 
explore desirable properties of the dynamics in each isolated element. The 
next Section 13.31 then uses centre manifold theory, based upon these results, 
to support the modelling of fully coupled dynamics. 





A 
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Theorem 2 (orthogonal eigenmodes) Consider the operator L on the 
right-hand side of the PDE with boundary conditions fllOI) -( TTT]) .- eigen- 
modes corresponding to distinct eigenvalues are orthogonal. 

Proof: Adapt the classical proof in many undergraduate texts. Let 
u and V be eigenmodes of C in the jth element corresponding to eigenvalues 
A and [J., respectively. For slight simplicity use f to denote fj. Since £u — Au 
and £v = (J.V , 

= (v, £u — Au) — (u, Cv — (xv) 
(integrating by parts) 

= [fvux - fuvj + [fvux - fuvj + ( (J, - A) (u, v) 

= f (Xhi )v(Xj+i )u,(Xj+i ) - f (Xj+i )u(Xhi )Vx(Xj+i ) 

- f (X+)v(X,+)u,(X,+) + f (X+)u(X+)v,(X+) 
+ f (Xr)v(Xr)u,(Xr) - f (Xr)u(XrK(Xr) 

- f (Xj_i )v(Xj_i )uJXj_i ) + f (Xj_i )u(Xj_i )vJXj_i ) 
+ (M-- A) (u,v) 

(by ([TT]) . u(Xj'^) is continuous and u(X-j-i-i) = u(Xj), 

and similarly for v) 
= v(Xj)[f(Xj+i)uJXj+i] -f(X+)u,(X+) 

+ f(Xr)u,(Xr)-f(Xj_i)Ux(Xj_i)] 
+ u(Xj) [ - f (Xj+i K(Xhi ) + f (X,+K(X+) 

-f(XrK(Xr)+f(Xj_iK(Xj_i)] 
+ (M-- A) (u,v) 

(by the boundary condition ffTOl) on gradients) 
= (M-- A) (u,v) . 

Hence, for distinct eigenvalues, A 7^ p., the corresponding eigenmodes must 
be orthogonal, (u, v) = . I|k 

Theorem 3 (reality) Consider the operator L in PDE ([9]) wzi/i boundary 
conditions (ITOl) -(|TTi) .• zfo eigenvalues are all real. 

Proof: Again adapt the classic proof. Let u denote any eigenmode 
corresponding to any eigenvalue A on the jth element; they are potentially 
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complex. Let v and (x denote the complex conjugate of u and A, respectively. 
Since the pde and boundary conditions fllOp - fllip have all real coefficients, 
these complex conjugates must also be an eigenmode /eigenvalue pair. The 
derivation within the proof of Theorem [2] then establishes that 

iyi-A) (u,v)=0. 

Here (u, v) = |up dx which is necessarily non-zero and hence |j. = A . But 
10. and A are complex conjugates, so any eigenvalue A must be real. 4 



Theorem 4 (spectral gap) Consider the operator C in PDE (Q with bound- 
ary conditions (fTOj) - (fTTj) .• when the coefficient function fj(x) is bounded above 
zero, fj,min — KiiUxeE- fj(^) > 0, then there is a zero eigenvalue and all other 
eigenvalues are negative and bounded away from zero by an amount propor- 
tional to f j,min/(Xj+1 ~ 

Proof: Let u denote any eigenmode corresponding to any eigenvalue A 
on the jth element. Thus Au = [f juJx • Multiply this ODE by u and integrate: 



A 



u^dx 



u[fjuJxdx 



f jU^ dx 



= fj(Xj+i)u(X^+i)u.(Xj+i) -fj(X,+)u(X,+]uJXt) 
+ f j (Xr)u(Xr)u,(Xr) - f ^ (Xj_i )u[X^^^ ]u,(Xj_i ; 

f dx 



(by ( |TTi) . u[Xf] is continuous and u(Xj-ti) = u(Xj)) 



= u(Xj)[f^(Xj+i)Ux(Xj+i) -fj(X^)uJX^) 

+ fj(Xr)u,(Xr) -fj(Xj_i)u,(Xj_0] - 
(by the boundary condition (ITOl) on gradients) 
f jU^ dx . 



f jU^ dx 
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Hence, for coefficient functions fj bounded above zero, the right-hand side 
is non-positive and hence so must all the eigenvalues A. The eigenvalue zero 
corresponds only to solutions that are constant on the jth element, — 

Now prove that all other eigenvalues are bounded away from zero by 
using a bound from the constant coefficient case of Section 13.11 Constrain 
the magnitude of the eigenmodes by u^dx = 1 . Then use the above 
identity to bound the eigenvalue 



A = 



f dx < — f j 



) ,mm 



u^dx. (14) 



Thus minimise u^dx subject to u^dx = 1 and u dx = (as other 
eigenmodes are necessarily orthogonal to the constant eigenmode) and the 
boundary conditions fllOp - (ITT]) . Using Lagrange multipliers [i and ■y, stan- 
dard Calculus of Variations asserts this minimum occurs for functions u(x) 
satisfying the Euler-Lagrange equation "v + 2(j.u — 2uxx = . For simplicity, 
let £, = X — Xj and Xj-^i — Xj = ±h, + h. Three cases arise: 

u, = -Mc^>0 =^ u= Acoshk£,-FBsinhk£,-F Csinhk|£,|; 

^ = ^ u= + A + B£, + C1^!; 
M, = -k^<0 =^ u = — Acosk£, + Bsink£,-F Csink|£,| . 

In the ffist case (hyperbolic), substituting u into the boundary conditions (fTOll - 
( ITT|) and the orthogonality condition gives four linear equations for A, B, C 
and "V which have nontrivial values only when the determinant 

sinh kh \ cosh kh — cosh khl = . 

k 

This has no real solutions for k. In the second case, substituting leads to the 
determinant 4h,^(h,+h,)(h,— H) , which also cannot be zero for non-degenerate 
grids. In the third case (trigonometric). A, B, C and y have nontrivial values 
only when the determinant 

— sin kh r cos kh. — cos khl = . 
k 

This only has solutions for finite k — for non-degenerate grids the smallest 
k = Tc/h — hence u^dx is bounded away from zero by an amount propor- 
tional to h^^ oc (Xj+i — Xj_i )^^, and thus so are the negative eigenvalues. ^ 
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Consequently the spectrum of the hnearised dynamics (pi)- (fTTi) on the 
jth element is qualitatively like that of the constant coefficient case (1131) . In 
particular, there always exists a spectral gap between the zero and the other 
eigenvalues. 

3.3 A relevant slow manifold discretisation exists 

Recall we aim to rigorously support discretisation of the nonlinear dynam- 
ics of the self-adjoint nonlinear general reaction diffusion equation ([7]). For 
definiteness we seek spatially periodic solutions, u(x + L,t] = u(x,t) for 
some period L, and divide the domain into m overlapping elements as shown 
schematically in Figure [T] (with X^j^ — Xq = L). The subgrid fields Uj(x, t) in 
each element satisfy the pde ([8]) and are coupled by the conditions ([5])-([6]). 
Then Theorem H] proves that linearised about any of the piecewise constant 
solutions of the subspace Eq there is a spectral gap in the dynamics of the 
reaction diffusion pde (^^. Consequently centre manifold theory [SI [221 ^-g-] 
asserts the following. 

Corollary 5 (slow manifold) For sufficiently smooth reaction g and diffu- 
sivity f in some neighbourhood of the subset o/Eq for which fj^min are bounded 
above zero: 

1. there exists a (m + 2) dimensional slow manifold Aio of the subgrid 
PDE ((Hj) coupled by ([5|)-([6|), one dimension for each element, and one 
dimension each for parameters y and a; 

2. the slow manifold Aio may be parametrised by any reasonable mea- 
sure Uj of the field in each element, that is, the slow manifold and the 
evolution thereon may be written, for some Uj and g^, ] = ],... ,vrL , as 

Ui = Uj(U, X, a,y] such that Uj = -^-^ = gj(ll, cx,y) ; (15) 

dt 

3. the dynamics on Mq is 'asymptotically complete' fSW in that for all 
solutions of the subgrid pde ([Hj) coupled by ([S|)-([SD from initial con- 
ditions in some neighbourhood of Aio, there exists a solution of the 
slow manifold model f|T5|) that is approached exponentially quickly in 
time — roughly at a rate minj{fj i„in/(Xj_|_i — Xj_T)^}; 
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4- the order of error of an approximation to the slow manifold }Aq and its 
evolution, fll5l) . is the same as the order of the residuals of the governing 
PDE dH]) and coupling ([5]) -([6]) when evaluated at the approximation. 

The evolution (fT5l) on the slow manifold M.o forms the discrete model 
of the dynamics of the PDE ([8]) coupled by ([5])-([6]). In principle, such a 
model is an exact closure for the discretisation in that the model tracks the 
evolution from general initial conditions [Corollary . However, we hardly 
ever can construct an exact slow manifold. Nonetheless, computer algebra 
readily constructs the slow manifold and its evolution to a controllable order 
of accuracy [Corollary [51H] . Then evaluating the model for full coupling, 
y = 1 , generates a model for the dynamics of the physical reaction diffusion 
PDE ([7]). Section provides evidence that y = 1 is within the finite domain 
of validity of the slow manifold Mq. 



4 Application to discretising diffusion 

This section briefly describes three interesting applications of the preceding 
Corollary [5l Section 14.11 investigate the simplest case of linear diffusion in 
order to show how the slow manifold varies with coupling parameter y. The 
interest is to see how the approach generates classic interpolation for the 
subgrid fields and classic finite difference rules for the evolution. Section 14.21 
explores nonlinear diffusion to demonstrate that this approach supports spe- 
cific discretisations for nonlinear problems and that, through using the cou- 
pling conditions ([5])-([6]), the discretisation preserves the self-adjointness of 
the original physical system. Lastly, Section 14.31 illustrates one way to ac- 
count for domain boundary conditions other than periodic, and also verifies 
convergence in the coupling y in one particular example. 



4.1 Linear diffusion 

The simplest application of the Slow Manifold Corollary [S] is to forming 
discrete models of linear homogeneous diffusion as governed by the PDE 

As described in Section |3l embed the L-periodic spatial domain into dynamics 
on m overlapping elements. Corollary [S] guarantees there exists a relevant 
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discrete, slow manifold, model of the diffusion fll6p . 

Computer algebra [Ml §2] readily constructs approximations to the slow 
manifold; one may check the approximation by confirming f|T7|l - f|T8|) satisfies 
the governing equations, f|T6|l and ([5])-([6]), to the specified order of error. Let 
the grid be uniform, AXj = h, and define the grid values Uj(t) = Uj(Xj,t) . 
Computer algebra finds the subgrid, intraelement field is 

Uj = [1 +y(£,^5 + 11^16^) +7^(1^25^- I|£,|52]]Uj + 0(y3) . (17) 

in terms of the subgrid variable £, = (x — Xj ) /h, , and the centred mean \x 
and difference 6 operators, |J.Uj = (Uj+1/2 + Uj_i/2)/2 and 6Uj = — 
Uj-i /2 • Reassuringly, when evaluated at full coupling y — 1 the terms linear 
and quadratic in the coupling parameter y form classic linear and quadratic 
interpolation, respectively, between the grid values Uj. The corresponding 
evolution on the slow manifold is the discretisation 



21 - 525y + 448y2 - 1 SOy^ 



1680h2 



V6^Uj + C(y^6^°). (18) 



Evaluated at full physical coupling y = 1 , this model recovers the classic 
centred finite difference formula for the discretisation. Computing to higher 
orders in coupling y, gives more and more terms in the classic formula. 

This approach recovers classic formula in such simple linear dynamics. 
However, the theoretical support is different: centre manifold theory applies 
at finite element size, to guarantee a relevant model for all initial conditions 
in some finite domain, and, after initial transients decay, for all times. The 
only approximation is the error incurred by the truncation of the description 
of the slow manifold in coupling parameter y. 



4.2 Nonlinear diffusion 

Nonlinear diffusion has many applications and is of continuing interest [21 
SSI e.g.] As a specific example application, consider nonlinear diffusion gov- 
erned by the pde 

9u 9 / 9u\ 



9t 9x1"^ 9^ J- ^''^ 
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As described in Section [31 embed the L-periodic domain into dynamics on 
m overlapping elements, with the coupling conditions on the gradient im- 
plemented with nonlinear diffusivity f(x, u, Ux) = u. Corollary [5] guarantees 
there exists a relevant discrete, slow manifold, model of the nonlinear diffu- 
sion (fT9D. 

The rate of attraction to the slow manifold is no longer a constant: instead 
it is now proportional to some measure of the smallest amplitude of the 
initial field. Consequently, high order asymptotic approximations to the slow 
manifold are replete with divisions by grid values Uj. Such divisors reflect 
the nonlinear diffusion and suggest that the domain of attraction of the slow 
manifold model reduces for fields u of small magnitude. 

Computer algebra [511 §3] readily constructs the subgrid field of the slow 
manifold: 

Uj = [^+y£,^5 + ly(^-y)|£,|62]+y2^£,262]Uj 

+ 0(y^). (20) 

The first line is identical to that for linear diffusion, (IT7|) : thus the second 
line is due to the nonlinearity in the diffusion. Unlike usual finite differences 
or finite elements methodology which imposes an interpolation between grid 
values, part of the value of this 'holistic discretisation' [33l e.g.] is that 
the subgrid field is constructed to satisfy the pde f|T9l) and thus generates 
accurate closures for the discrete model. The evolution on the slow manifold 
gives the discrete model 

Truncated to errors O(y^) then evaluated at full coupling, y = 1 , this is 
the introductory model ([1]). To errors of order O(y^), the discretisation (I2T!) 
of the nonlinear diffusion is the same as the discretisation fll8p of linear 
diffusion, but applied to -^U? instead of to Uj. This is reasonable since this 
continuum nonlinear diffusion operator (uUx)x = (^"U-^jxx- This nontrivial 
correspondence, not imposed at all but instead a natural closure from the 
intricate subgrid scale interactions, confirms this approach to discretisation 
is sound. 

However, such a very close correspondence breaks down at O(y^) when 
divisions by \lj±] start invading the slow evolution fl^Tl) of the nonlinear 
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diffusion. Subtleties in such higher order subgrid scale interactions result in 
such more complicated discretisations. 

Nonetheless, the equivalent differential equation of fl2T|) is 

— - y^A fu— ^ +hV(i -y)^"^^— (^3— ^ + 

9t dx \ dx J 12 9x \ 9x 9x^ 9x^ y 

+ 0{}x\y'), (22) 

which is in conservative form no matter what order we truncate the anal- 
ysis in coupling parameter y. The coupling conditions ©-(jo]) do preserve 
conservation. 

4.3 Dirichlet boundaries on one element 

In every other section we explore dynamics far away from physical boundaries 
by assuming spatial periodicity. Conversely, this section explores the extreme 
case of precisely one element between physical boundaries forming that one 
element, say the element is non-dimenisonalised to —1 < x < 1 . This ex- 
treme case empowers us to compute to high order and show convergence in 
the 'coupling' parameter y, as well as illustrating the ease of incorporating 
physical boundary conditions on the global domain rather than assuming 
periodic conditions for the m elements as used elsewhere. 

For definiteness, suppose the physical boundary conditions for the exam- 
ple PDE of nonlinear diffusion (IT^ are the Dirichlet conditions that u = 
at X = ±1 . Implement these Dirichlet conditions on the one element [—1 , 1] 
via the adaptation of the coupling conditions ([5])-([6]) to 

u(0+, t]ux(0+, t) - u(0", t]ux(0", t) 

= Y [u(1 , t]ux(1 , t] - u(-1 , t]ux(-1 , t]] , (23) 
u(±1,t) =y'u(0,t) and u(0", t) = u(0+, t) . (24) 

When coupling parameter y = (y' = 1) the element is isolated from 
the physical boundaries and the spectrum of the dynamics on the one el- 
ement are as for the previous Section 14. 2[ Thus there exists a slow manifold 
parametrised by y and the mid-element value Uo(t) = u(0, t) : the slow man- 
ifold may be described by u(x, t] = Uo(x, Uo,y) such that Uq ~ g(Uo,y) . It 
is exponentially quickly attractive, roughly at a rate oc Uq . 
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Straightforward adaptions of the computer algebra [M] construct the slow 
manifold and the evolution thereon to high order in the coupling parameter y. 
The slow manifold is 

Uo = Uo[l -yix|+y2(lx|-x2)+y3_^(|x|+9x2-10|xp)+O(y4)] . (25) 

Observe the y^ and y^ subgrid structures give classic linear and quadratic 
interpolation when evaluated at the physical y = 1 . It is the O(y^) terms 
that begins to account for the nonlinear nature of the diffusion, hence pro- 
vide correct subgrid structures, and consequently provide a sound closure for 
the macroscale, one dimensional model on the element. The corresponding 
evolution on the slow manifold is 

'-^0--Uo[y +27 ~487 "96^ + 8647 + 8647 + (,7 jj • 

(26) 

This predicts that from all nearby initial conditions, and apart from expo- 
nentially quick transients, solutions decay algebraically in proportion to 1 /t 
as 'material' u diffuses through the physical boundaries at x = ±1 where the 
diffusivity is zero but the flux is not. 

With just one 'amplitude' Uo computer algebra generates approxima- 
tions to (9(y^°) within a minute. Then generalised Domb-Sykes plots [26l 
Appendix] estimate the location of the convergence limiting singularity. Fig- 
ure [2] strongly suggests convergence in coupling parameter y with radius of 
convergence 1.18 due to a complex conjugate pair of (logarithmic) singulari- 
ties in complex y at an angle of 54° to the real y axis. These plots provide 
good evidence that evaluation at y = 1 of the power series' is convergent for 
at least this simple case. 

5 Disjoint lattice elements do not couple 

Research into spatio-temporal dynamics commonly invokes a discrete lat- 
tice [HI [ini EH IMl US, e.g.]. Section O introduces rigorous support for a 
transformation from fine scale lattice dynamics to a coarser scale lattice dy- 
namics. Similar to the continuum case of Sections [SHU the transformation 
requires embedding the dynamics in a system of twice the dimensionality, 
before reducing the dimension by a factor of four, to achieve an overall halv- 
ing of the dimensionality. The resultant model then resolves dynamics on a 
lattice with twice the grid spacing. This section argues that it is difficult. 
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(b) l/n^ 

Figure 2: power series to errors (9(y^°) generate these generalised Domb- 
Sykes plots [261 Appendix] that, when extrapolated to l/n = , strongly sug- 
gest the power series in coupling parameter y has radius of convergence 1.18. 
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G • • © 

Uo Ui U2 U3 U4 Us 

Figure 3: rewrite the dynamics of u as the dynamics of v on two elements by 
renaming variables: the solid discs correspond to differential equations, and 
the open circles correspond to algebraic coupling equations. 

if not impossible, to achieve such net halving of the dimensionality without 
the initial doubling via the embedding. 

As the most basic, but key, example of self-adjoint dynamics on a lattice, 
consider the discrete diffusion equation nondimensionalised to 

Ui = Ui_i - 2ui + U|+i , (27) 

on equi-spaced grid points Xi = ih, . This section seeks to construct a sound 
coarse scale model for these dynamics, but fails because of interesting reasons 
that inspire the next sections. 

For almost extreme simplicity, suppose the discrete diffusion fl271) applies 
on a small finite domain at just four lattice points, i — 1,2,3,4, with insu- 
lating Neumann-like boundary conditions provided by 

Ui — Uo = U4 — U5 = . (28) 

We seek a model in just two dynamical variables, for this system with four 
dynamical variables. Because the dynamics are so low dimensional, and 
linear, we reasonably explore all options. 

Divide the domain into two elements, the first containing ui and U2 and 
the second containing U3 and U4. As shown in Figure [3l to connect the 
two elements introduce new names V3 = V5 = U3 and V2 = V4 = U2 for 
the middle two dynamical variables. For convenience, rename the others 
variables Vq = Uq , Vi = ui , Vg = U4 and V/ = U5 . Then write the discrete 
diffusion fl27|l and insulating boundary conditions fl28l) . together with the 
interelement coupling identities, as the differential-algebraic system 

Dv = Lv , (29) 
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where v = (vq, . . . , V/], D = diag(0, 1 , 1 , 0, 0, 1 , 1 , 0) and the matrix 



- 1 ] 
1-21000 

1-2 1 

1 0-10 

~ -1 ^ 

1 -2 

^ 

_ 

where I explain the import of the itahc entries shortly. The spectrum of the 
linear discrete diffusion fl29|) . from det(L — AD), is {0,—2 + V2,—2,—2 — \/2.}. 
We construct a long term model from the two slow modes corresponding to 
the two eigenvalues nearest zero. Thus we seek a model with just two dy- 
namical variables that systematically track the amplitude of the two slowest 
modes. 

Centre manifold theory provides rigorous support for such low dimen- 
sional modelling [8], [22], e.g.]. But, analogously to the continuum analysis of 
Sections [SHU we need to artificially modify the linear operator L to have two 
eigenvalues of zero, then implement a homotopy (a smooth path) to recover L 
and the original dynamics [33l [35l e.g.]. I argue this is impossible, and hence 
we need the more complicated embedding of the next section. 

To see what freedom we have available, first identify the aspects of L that 
cannot be changed. We posit that the evolution equations, corresponding to 
the second, third, sixth and seventh lines of L cannot be changed as they are 
to encode the microscale dynamics (127|) : if we need to modify the microscale 
dynamics, then any embeddings we find are almost certainly problem specific 
and thus not of general power. To maintain the self-adjoint symmetry of L, 
this then fixes Li_2; ^4,3, ^5,6 and Lg j to be one, and many other entries to 
be zero. The entries to vary are those in italics in (150]) . Second, of these, set 
1-1,8 — hs,! = to avoid excessively nonlocal equations. Third, preserving 
the zero eigenvalue of Ut = constant, the conservation mode, we need each 
row sum to remain zero. Lastly, as well as self-adjoint symmetry, we require 
isotropy, left-right symmetry. These four requirements result in three degrees 
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(30) 
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of freedom spanned by the three matrices 

0" 







0' 





0_ 

' 







1 0-1' 



-10 1 _ 

1 " 


1 
-1' 


-1 _ 

Thus the (nearly) most general self-adjoint, isotropic, conservative, matrix of 
the discrete diffusion dynamics is L-|-'yiYi+-y2Y2+'y3Y3 for some constants y i , 
yz and 1^3. Its spectrum is given by the zeros of the characteristic polynomial 

det(L + yiYi + 'y2Y2 + ysYa - AD) = CiA + CaA^ + CsA^ + C4A4 , (31) 

for some coefficients depending upon y^,y2 and yy. for example, 

C4 = -(1 -2ij2)(1 +4ij^ + 4y2y3-2yi +2yi-y2). 

Recall that we seek to create a homotopy from dynamics with a double zero 
eigenvalue on two decoupled elements, iji = 1 and 1^3 = (the 'green line' in 
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-0.2 

-0.6 

-1.0 



-1.0 



Figure 4: red: iso-surface of 'infinite' eigenvalues from the coefficient C4 = 
of the characteristic polynomial (l3Til . Blue blob: the origin indicating 
the original diffusion dynamics (I3UI) . Green line: two decoupled elements. 
There appears no route from green to blue without encountering 'infinite' 
eigenvalues, red. 
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Figure H]), to the original dynamics with 1^1=^2=^3 = (the 'blue blob' 
in Figure H]). Figure H] indicates that it is impossible to create a homotopy 
from the two decoupled elements to the original dynamics without crossing 
the red surfaces. The red surfaces in Figure H] are the surfaces of C4 = 
in the charateristic polynomial fl3Tl) . and hence represent neighbourhoods 
where eigenvalues become infinitely large, both positive and negative. We 
cannot create a smooth homotopy from two decoupled elements to the cou- 
pled original dynamics through such neighbourhoods. This failure with just 
two elements suggests that, within the class of self-adjoint, isotropic, diffusion 
dynamics, we cannot artificially divide a domain into disjoint elements. 

Thus the next section proceeds to explore embedding lattice dynamics 
onto overlapping elements analogous to the overlapping elements used for the 
continuum dynamics of Sections [2H1] and for other multiscale approaches [131 
SQlIISl e.g.]. 

6 Overlapping elements preserve self-adjoint 
dynamics 

This section explores how to transform the discrete dynamics of variables Ui(t) 
on a fine grid of spacing h, into discrete dynamics of variables Uj(t] on a 
coarser grid of spacing H = 2h. 

Figure [5] schematically shows that the jth element stretches from a neigh- 
bourhood of Xj_i to a neighbourhood of Xj+i . Figure [5] also shows each 
element is divided into two halves, and the variables duplicated so that, no- 
tionally, Uij = Vj^i^g = Vj^j = V],4 = Vj+1,0 and U2j+i = Vj_i,7 = Vj3 = Vj,5 = 
Vj+i 1 . Embed the fine grid dynamics in these overlapping elements in a space 
of double the dimensionality by treating Vj_i _6 and Vj 2, and Vj ^ and Vj+i j as 
independently evolving variables. 

In the embedding space we use the inner product (v, w) = J^ - -i^Vj ^wj i. 

6.1 Full coupling rules 

Analogous to the continuum dynamics of Sections [2H11 we need to find rules 
to couple neighbouring overlapping elements that are sufficiently generic that 
they are useful for modelling a wide variety of self-adjoint lattice dynamics. 
Here we focus on the basic case of the coarse grid modelling of the discrete 
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Vj-1,0 "^3-1,3 
G • • © 

■^j-I.S "^3-1,7 

G ■ • O 

Vj,0 Vj,^ Vj_2 Vj_3 
G • • © 



(j — 1 )th element 



G • • Q 



jth element | Vj,4 Vj,5 Vj^g 

(j + 1 ]th element 



Vj+1,0 "^3+1,2 Vj+l_3 

G • ■ © 

Vj+1,4 ■^3+1,5 "^3+1,6 Vj+l_7 
G ■ • 

U2j-4 U2j-I •^Ij-f \!;23-f \!;23-l ^2j+^ ^2j+^ ^Zj+f \!;2j+I ^Ij+f 



Figure 5: to transform dynamics from the fine grid, bottom, to the coarse 
grid, top, rewrite the dynamics of Ui as the dynamics of Vj ^ on overlapping 
elements, here see three consecutive elements (blue, black and magenta), by 
duplicating and renaming variables: the solid discs correspond to differential 
equations showing the number of dynamic variables is doubled, and the circles 
correspond to algebraic coupling equations. Each element has two halves, 
also shown separated for clarity. 



diffusion dynamics (127|) : equation (!27|) applies at the internal dynamic vari- 
ables Vj 1, Vj_2, Vj 5 and Vj g of each element, the discs in Figure [51 Now take 
up the challenge of using the variables Vj,o, Vj,4, Vj,5 and v^j, the open circles 
in Figure [5], to couple together not only the elements but also the two halves 
of each element . 

Again analogous to the continuum dynamics of Sections [2HU where rig- 
orous theorems are invoked we require that the lattice dynamics are periodic 
in the fine grid with period 2m in i, and thus periodic on the coarse grid 
with period m in j. 

Let us identify the possible domain of evolution and coupling rules. Fol- 
lowing Section [5], one evolution rule within each element is 

DVj = LVj (32) 

where Vj = [vj^, . . . ,Vjj], D — diag(0, 1 , 1 , 0, 0, 1 , 1 , 0) and equation (!30!) 
gives the matrix L. As in Section [5l to preserve the self-adjoint isotropic 
conservation of the dynamics the upright entries in the matrix L are fixed; we 
can only modify the italic entries. As in Section one set of possible changes 
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to L is spanned by Yi, Y2 and Y3. This set changes only internal interactions. 
There are also three basic matrices that couple elements together preserving 
self-adjoint isotropic conservation, with the additional constraint that the 
coupling has to be 'local' on the grid. Let e± denote the shift operators from 
one element to its neighbours: define £±Vj,i = Then write three basis 

matrices for nearest neighbour coupling changes to L as 
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For example, the non-zero elements of Y4Vj are vj o ~ '^j-1,7 E^-iid Vj 7 — Vj+i _o 
which connects end values of neighbouring overlapping elements. Thus we 
explore the self-adjoint isotropic conservative dynamics of 

DVj= (l + Y_v,Y^v^, (33) 
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for some coefficients yi that we are free to choose. 

When fully coupled, the element dynamics (15^ must reduce to that of 
the discrete diffusion (1271) . These are most compactly expressed in terms 
of the fine grid shift operator e± defined as £±Vj i = Vj_i±i . (When used 
appropriately, the coarse grid shift e± — £^ .) Then the fine grid diffusion, 
expressed as Ut = (£+ — 2 + £_]ut, has 'eigenvalue' (£+ — 2 + £_) correspond- 



ing to 'eigenvector' u 



where all are interpreted 



in an operator sense. For the element dynamics (133|) to reproduce these dy- 
namics exactly, the dynamics must have the same 'eigenvalue' (£+ — 2 + £_) 
but corresponding to the 'eigenvector' = £^(1 , £+, £+, £^, £+, £+, £+, £+). 
Substituting these into equating coefficients of £± in all components 

gives a set of equations that uniquely determine y-\ =1^5 = ^6 = 1 and 
yz = Va = V4 — . Thus the unique operator that gives the correct evo- 
lution of the diffusion fl27|) with the correct subgrid microstructure for the 
diffusion is 
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(34) 



In a homotopy from a useful base for centre manifold theory, we must end the 
homotopy at this operator for the fully coupled dynamics on the elements. 



6.2 A homotopy connects elements 

Now seek a base of decoupled elements from which a slow manifold may be 
constructed to model the coarse grid dynamics. Perhaps the closest operator 
to the fully coupled Li of (IMl) is simply to change the inter-element coupling 
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shift operators £± to be ones: that is, define 



Lo = L + Yi+Y3- 
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(35) 



Then straightforward algebra derives that the spectrum of Dvj = LqVj is 
{0, —2/3, —2, —4} for each of the m decoupled elements. The zero eigenvalue 
with all the rest negative implies. Section 16.31 that there exists a relevant 
slow manifold model [H [221 e.g.] which we can construct with one dimension 
for each element — the dynamics on the m-dimensional slow manifold forms 
the coarse grid model. 

But is Lo a good choice? and, is it the only choice? The spectrum of the 
homotopy answers. Create a general homotopy from decoupled elements to 
fully coupled by defining the convex combination 

I-Y = (1 -y) [Lo + yiYi +y2Y2 + y3Y3] +yl-i , 

where parameter y morphs the operator from the decoupled case, y = , to 
the fully coupled case, y = 1 . The characteristic polynomial, det(LY — AD), 
is too hideous to record here, but computer algebra derives it easily. Then 
computer algebra iteration finds asymptotic approximations to the eigenvalue 
near zero: 

A = y'±{li _ 2 + 4) [1 + I(y 1 - 181,2 - lOija)] +0{y' + \y\^) . 

The factor j^{e^ — 2 + e^), although unexpectedly involving double shifts 
over the coarse grid, is exactly what we want as to leading order it is the 
centred second difference 6^ of the fine grid diffusion dynamics fl27|) . However, 
the term linear in yi ruins this identity and so we choose y-[ — ^8y2 + 
10ij3 . Similar computer algebra to higher order in |y | indicates we also need 
ys = —31,2, and analysis to higher order in coupling parameter y indicates 
that then we need 1,2 = . To best match the dynamics of the discrete 
diffusion (1271) we choose y-\ =1,2 = ^3 — 0- Consequently, accuracy of 
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the coarse grid model for just simple diffusion requires us to use the convex 
combination operator 



as the homotopy to smoothly connect the decoupled matrix fl35|) with the 
fully coupled operator flMl) . 

Inhomogeneity requires generalisation Linearisation in Section [6731 leads 
us to consider the class of inhomogeneous, discrete, reaction-diffusion equa- 
tions on the fine grid of 

Ui = 5(fi6ui) -hagiUi = ft_mi_i - (ft_^ +fi+^]ui-Ff^_^mi+i -agiUi , (37) 

for some spatially varying 'reactions' Qi and 'diffusivities' f^^i governing the 
dynamic exchange between Ui and Ui±i. Equation fl37|l describes relatively 
general self-adjoint dynamics on the fine grid. 

Embed the fine grid dynamics fl37j) into the dynamics on the finite ele- 
ments shown in Figure E] and generalising the self-adjoint, consistent opera- 
tor (1361] by considering the dynamics 



y'Lo + yLi 



where normally y ' = 1 — y , 



(36) 



(38) 



where the linear 'reaction 



9j = (0, 92^-1 Vj,i,g2j_^Vj,2, 



0,0, 92^+1 Vj,5,g2j+|Vj,6,0) , 



and where the diffusivity matrix 




for the four sub-blocks 
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Fi,2 = diag[(y'+y£_)f2j,0,0,f2j(y'+y£-)], 
F2,i = diag[f2j(y' + y£+],0,0,(y' + y£+)f2j], 



recalling that £±f2jVj^ — f2j±2Vj±i,i. The operator on the right-hand side 
of fl38|) is self-adjoint: the only subtlety arises via the operators coupling 
elements. For example, suppose the operator £±f2j occurs at element (i, k) 
in Fji then the inner product (w, Fv) has the components 



which suitably corresponds to the (k, i) elements of Fj in the inner prod- 
uct (Fw,v). Thus the dynamics of the embedded, coupled, finite element, 
system (1381) preserves the self-adjointness in the fine grid inhomogeneous 
dynamics flTTI) in a manner consistent with the required homotopy fl5B]) of 
homogeneous dynamics. 

6.3 Centre manifold theory supports a coarse grid model 

We have arrived at a separation of the self-adjoint, isotropic, diffusion dy- 
namics 0271) into elements, drawn schematically in Figure El with coupling 
between the elements that ranges from decoupled to fully coupled, and pre- 
serves the self-adjoint nature of the linear dynamics. We now invoke centre 
manifold theory [H [22], e.g.] to rigorously support coarse grid modelling of 
nonlinear lattice dynamics. The nonlinear dynamics are introduced, embed- 
ded into elements, and then linearised to connect to the preceding results. 

In analogy with the continuum dynamics of ([7]), here we consider the class 
of fine grid dynamics expressible by the general nonlinear, local interaction, 
lattice rule 



for sufficiently smooth 'diffusivities' f and 'reactions' g. For definite theoret- 
ical statements, suppose the fine grid lattice and the dynamics fl40|] on it are 
2ra-periodic in i. 






6[f(i, [lUi, 5ui)5ui] + ag(i,Ui, |j.6ui) 



(40) 
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Embed the fine grid, 2m-periodic, nonhnear reaction-diffusion dynam- 
ics (1401) into the higher dimensional dynamics of Vj within the corresponding 
m elements of Figure [5] by a nonlinear version of fl38|) . To help write the map 
from the fine grid to the overlapping elements let i' — [i — j) — sign(i — |), 
then the nonlinear fine grid dynamics (HOl) in each element are 

Vj,i = 6[f2j+i'Svj,i] + <^g2j+i', 1=1,2,5,6, (41) 

where fzj+i' = f (2j (xvj^i, 5vj_t) and gij+i' = g(2j Vj_i, (J.6vj,i). Control 
the coupling of these to neighbouring elements by the conditions 

Wi^^hi + iY + y£±) [f2j6v. z] = , i = { , (42) 

f2j(y' + y£+]vj,o + izMii - f2j(y' + yi-]v^,7 = , i = f, | . (43) 

The embedded element dynamics (14T]) - (143|) then has a useful m dimensional 
subspace Eq of equilibria: a = y = and Vj ^ = Uj constant in each of the 
m elements. 

Linearise the dynamics about each equilibria in Eq to obtain (1381) without 
reaction, oc — , and with diffusivites f2j+i = f (2j U,-, 0). As the coupling 
parameter y — , each element is isolated. When the diffusivities ft are 
constant, each element has spectrum proportional to {0, —2/3, —2, —4} — the 
zero eigenvalue corresponds to Vj t being constant in each element. Elemen- 
tary algebra shows that provided for all j 

fj>0 and 2f2j(f2^-2 + f2j+2)-f2j-2f2j+2>0, (44) 

then the spectrum remains as one zero eigenvalue with the other three being 
negative. Consequently, centre manifold theory assures us of the following 
corollary P [221 e.g.]. 

Corollary 6 (slow manifold) Provided (H4l) . in some finite neighbourhood 
of the subspace Kq: 

1. there exists a (m -|- 2) dimensional slow manifold M.q of the nonlin- 
ear, fine grid, element dynamics (14T]) - (l43l) . one dimension for each 
element, one for the inter-element coupling parameter y, and one for 
the amplitude parameter a of the reaction; 
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2. the slow manifold may be parametrised by any reasonable mea- 
sure Uj of Vj i in each element, that is, the slow manifold and the evo- 
lution thereon may be written, for some Vj and gj, j = 1 , . . . , m , as 

Vj=Vj(U,y,a) such that Uj = = gj(U, y, a) ; (45) 

3. the dynamics on is 'asymptotically complete' in that from all 
initial conditions in some neighbourhood of Aio, there exists a solu- 
tion of ( H5|) approached exponentially quickly in time by the solution 

of m-m; 

4- the order of error of an approximation to the slow manifold Aio and 
its evolution, (1451) . is the same as the order of the residuals of the gov- 
erning dynamics ( 1411) - (H3l) . in the coupling parameter y and reaction 
parameter cc, when evaluated at the approximation. 



Basic example of linear diffusion The simplest example of the class 
of lattice dynamics to which Corollary [6] applies is the linear discrete dif- 
fusion (l27j) . Computer algebra [SU §4] readily iterates to asymptotically 
approximate the slow manifold model. Here we choose to parametrise the 
slow manifold of coarse scale dynamics in terms of the coarse variables 

Uj = ^(Vj,2+Vj3+Vj,4 + Vj,5), (46) 

which Figure \5\ shows to be estimates of the mid-element values of the fine 
grid variables. Executing the computer algebra deduces that the fine grid, 
intraelement structure is 



(1,1,1,1,1,1,1,1) 
+ |(-5,-3,-1, 1,-1, 1,3,5)^x6 

+ ^(3,1, 0,0,0,0,1, 3)(5^+l64) 



+ 1^(13,7,1, -5,5,-1, -7, -13)H63jUj + 0(y3), 

in terms of the centred mean p. and difference 5 operators on the coarse grid. 
The first line gives the piecewise constant basis for the slow subspace Eq 
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of equilibria. The second line gives a linear variation between neighbouring 
elements, the third line quadratic, and so on. Corollary E] assures us that no 
matter what the initial conditions for the fine grid diffusion fl27j) . exponen- 
tially quickly the system will settle onto the slow manifold with the above 
structure on the fine grid. 

The corresponding evolution of the coarse grid variables Uj is 



+ y'{¥' + ^s'^' + m,'^')h + 0{'Y'). (47) 



Obtain the lowest accuracy model from the first line, by neglecting terms 
C(y^), and then evaluating at the physical fully coupled case y = 1 : the 
model is U, = ^(Uj_2 ~ 2Uj + Uj+2) which is appropriate although sur- 
prisingly only involves every second coarse grid value. This 'surprise' was 
forecast in Section [6l2] by the leading operator eigenvalue of the operator Ly. 
Higher orders in coupling parameter y give coarse grid models with a wider 
stencil, and of more accuracy. For example, we see the accuracy through the 
equivalent fine grid expression of the coarse grid model (H7I) obtained by the 
operator identity that 5^ = 48^ + 6"* : in terms of fine grid differences 6, the 
coarse model fHTj) is 

= [yV + |y2(1 -y)(1 -3y]64] U^ + 0{&',y^) . 

This equivalent fine grid expression demonstrates that when evaluated for 
full coupling, y = 1 , the fourth order differences vanish to leave the correct 
equivalent model Uj ^ 5^Uj . Computer algebra [Ml §4] to high order in 
coupling y confirms that higher order differences in the equivalent fine grid 
expression similarly vanish. Thus the coarse grid model (147|) is an accurate 
closure for the coarse scales of the fine grid dynamics. 

As introduced in earlier work with non-self-adjoint coupling [35], such 
a mapping of lattice dynamics from fine grid scale to the coarser scale can 
be iterated and renormalised to cover step-by-step the wide range of length 
and time scales on a multigrid [3, HH e.g.]. Such step-by-step dynamical 
transformations could empower us, in future research, to carefully explore 
the development of emergent phenomena in nonlinear and stochastic systems 
over multigrids. 
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7 Three further apphcations indicates range 



The previous section concluded by modelling the dynamics of discrete dif- 
fusion. By itself, discrete diffusion is well understood. The value of the 
preceding section is that it empowers us to model more complicated dynam- 
ics from the same base with the same powerful centre manifold support. This 
section introduces three interesting applications. 

7.1 Reaction-diffusion lattice dynamics 

Centre manifold theory, as recorded in Corollary [6|, applies to nonlinear dy- 
namics such as the coarse grid modelling of the discrete reaction-diffusion 
equation [l6l e.g.] 

Ui = 6^Ui + a(ui - u?) , (48) 

where a(Ui — u?) is some example nonlinear reaction. 

Simple modifications of earlier code [SH §5] gives computer algebra that 
constructs the slow manifold model of the coarse grid evolution. Executing 
the code derives, for example, the coarse grid model 

ilj = |^(Uj+2-2Uj + Uj_2) + a(Uj-Uf) 

2 

+ ^(l4Uf - 10UjUf+i - 5Uf+T -4UjUH2 + 10Uj+iUj+2 

- 3Uf^2 - 10UjU^_i + 10Uj_iUj+i - 5Uf_, - 4U^Uj_2 

+ 10Uj_iUj_2-3Uf_2) +0(y',a^). (49) 

The first line, at full coupling y = 1 , is a classic model of the reaction- 
diffusion (148|) . The second and subsequent lines account for sub-element 
reaction and diffusion, interacting together and with neighbouring elements. 
Such terms are required for an accurate closure of the nonlinear fine grid 
dynamics on the coarse grid. 

7.2 Homogenisation 

A critical issue in material science is the effective large scale properties of a 
composite material with significant microscopic structure. A canonical prob- 
lem is the effective large scale diffusion through a domain with microscopic 
variations in diffusion coefficient [H SOI [D e.g.]. Here we transform diffusion 
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on a fine grid, with fine grid variations in coefficient, into diffusion onto a 
coarser grid. Here we provide a new and powerful view of the classic result 
that, to leading order, the coarse grid diffusion is a local average of the fine 
grid diffusion. 

Consider the fine grid dynamics PUI) . without reaction (x — , where the 
diffusion coefficient governing flux between fine grid points x^ i and X|_|_i is 
f i = 1 + e Kt . The parameter e controls the overall size of the variations in 
the diffusivity. For simplicity in construction and interpretation I treat the 
variations in diffusivity as small; that is, we construct the slow manifold as a 
power series in e. Computer algebra then finds the slow manifold as a power 
series in the strength e. Straightforward modifications to earlier computer 
algebra [Ml §6] finds, for the example ([2]) mentioned in the Introduction, 

= y'i^{/CMUj-2 - (/Cm + /Cj+OUj + /Cj+iUj+i} + o(y') , (50) 

where the coarse grid effective diffusivity coefficients 

ICj = ^ +el [k2]-2 + K2j-1 + K2j+1 + K2j+2] 

+ ^^16 [ ~ ('<2j-2 + K2j-1 - K2j+1 - K2j+2]^ 

-2(K2j-2-K2j-l)'-2(K2H2-K2Hl)^] • 

Because the element coupling preserves self-adjoint symmetry we find that 
the coarse grid, slow manifold model fl50l) is indeed self-adjoint with these 
particular effective diffusivities. 

Even more beautiful is that the effective diffusivity on the coarse grid 
is local: the diffusivity /Cj±i governing the flux between coarse grid points 
Xj and Xj-i-2 depends only upon the fine scale diffusivities between Xj and Xj±2, 
namely between Kjj and K2j±4. This beautiful feature is not built into the 
approach, but appears as a natural consequence of this scheme to preserve 
self- adjoint properties. 

Zigzag microstructure The specific zigzag microstructure = (— 1)^ 
is straightforward to analyse to higher order [Ml §6], and the resultant 
macroscale model compact enough to record. With zigzag microscale dif- 
fusivity, the diffusivity is uniform on the coarse grid: writing the coarse scale 
evolution as Uj = VXlj , the coarse grid diffusion operator 
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(51) 



As it should, the 1 — factor shows that the coarse scale diffusion is de- 
pressed by the microstructure, and indeed vanishes for e = ±1 reflecting the 
vanishing of the microscale diffusivity between every second fine scale grid 
point that occurs at e = ±1 . 

7.3 Pattern evolution 

An outstanding issue in modelling is the direct discrete, macroscale modelling 
of pattern evolution [HI [3], e.g.]. Currently, in fields such as fluid convection 
and in some reaction- diffusion equations, we first model the evolution of rolls, 
spirals and spots by variants of the Ginzburg-Landau pde derived assuming 
amplitudes and phases of the structures vary slowly in space. Second, we 
then discretise the Ginzburg-Landau pde for numerical simulations of the 
modulation of the pattern over large scales. The challenge is to analyse the 
original system dynamics and directly generate such a macroscale discreti- 
sation in one step. My first attempt to do this showed potential [Hlj, but 
the coupling conditions used therein did not preserve the self-adjoint nature 
of the dynamics and so the discretisation unsatisfactorily did not preserve 
required symmetries. Now with self-adjoint coupling conditions we return to 
modelling pattern evolution. 

Linear dynamics On a lattice, the simplest microscale pattern is perhaps 
the long lasting, zigzag mode of the discrete equation 



Now adapt the self-adjoint coupling conditions and analysis of Section [6] to 
derive a model of these zigzag mode dynamics on the coarser grid. Embed 
the dynamics on the elements shown in Figure \5\ with the evolution rule 



U| = — 4|j.^Ui 



(52) 



(53) 
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where the 'zigzag' operator 
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Then based about the decoupled case of y = , centre manifold theory [H [221 
e.g.] similarly guarantees the existence and relevance of a slow manifold 
parametrised by the amplitude of the zigzag mode local to each element. 
Computer algebra [Ml §7], modified from earlier code via spatial patterns 
changed to (—1 Y to account for the zigzag neutral mode, constructs the slow 
manifold model. The corresponding evolution of the coarse grid order pa- 
rameters, the amplitudes Uj, is exactly ( l47l) discussed before. The difference 
is that here the equation governs the local amplitude of the zigzag mode, 
rather than the local mean field. 



Nonlinear amplitude modulation Even more interesting is nonlinear 
lattice dynamics which centre manifold theory also supports. Modify the 
discrete equation (152!) by a local quadratic nonlinearity to 

Ui = — 4|j.^Ui + auf = — Ui_i — 2ui — Ui+i + au? . (54) 

What is the corresponding coarser grid equation governing the local am- 
plitude of the zigzag mode? In particular, does the quadratic nonlinearity 
stabilise or destabilise the origin? 

Simple modifications to the computer algebra then derives the coarser 
scale, modulation equation 

Uj = + - aylCSm^]{b%) + |a^Uf + ^(y^ + a^) . (55) 

Evaluated at y = 1 this coarse grid model predicts that the quadratic nonlin- 
earity in flMl) is destabilising as it generates the cubic growth term -^a^U? . In 
other problems, the coarse grid model will similarly form a discrete version of 
the Ginzburg-Landau equation. One might derive such a Ginzburg-Landau 
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equation by traditional modulation theory, but here analysis of the fine scale 
grid dynamics generates the novel term ay;^(6p.Uj)(6^Uj) which accounts 
for interactions over the scale of a few grid points and thus forms a more 
comprehensive closure. 

This analysis empowers us to derive discrete models of nonlinear lattice 
pattern dynamics without having to assume the infinite scale separation re- 
quired by traditional modulation theory. 

8 Conclusion 

This article uses centre manifold theory to further develop a novel approach 
to high quality coarse scale discrete models of nonlinear spatiotemporal sys- 
tems. The method is to divide space into overlapping finite elements with 
specially crafted coupling conditions. The innovative coupling conditions of 
Section [2] not only engender centre manifold support, Section |3l and assure 
consistency for vanishing element size, but also preserve the self-adjoint sym- 
metry which is often so important in applications. 

The first half of this article extracts accurate finite scale discrete lattice 
models from analysis of the infinitesimal scale modelling of dissipative pdes. 
A companion problem is to extract an accurate coarse scale lattice model 
from a fine scale lattice model. The second half of this article does this for 
the possibly extreme case when the coarse scale lattice is just a factor of 
two coarser than the fine scale lattice, and thus connects this approach to 
multigrid methods [7]. The approach is again based upon dividing the fine 
lattice into overlapping elements with specially crafted coupling conditions. 
Section O To supplement an earlier approach [53] , here the coupling condi- 
tions preserve self-adjointness in the dynamics as seen in the three example 
applications of Section [71 

Importantly, centre manifold theory [8], [221 6.g.] supports modelling for 
finite spectral gap, not just the infinite spectral gap required by other meth- 
ods. Furthermore, the support applies, in principle, to a finite domain and 
for all time; the only approximation is in the approximation of the slow man- 
ifold. The approach developed here for deterministic dynamics in one spatial 
dimension should be extendable to both higher dimensions and stochastic 
dynamics. 
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